Analysis_01_prone_incidence_risk_reliability
Chad H Hochberg and Anna Barker
2025-08-22
#Specify Project Location and File Type
project_location <- '~/Library/CloudStorage/OneDrive-JohnsHopkins/Research Projects/CLIF/CLIF_Projects/ProneProjects/ProningEpi/ProneSevereARF_Output'
file_type <- 'csv'
#Create Sub Folders within Project Folder
# Check if the output directory exists; if not, create it
if (!dir.exists(paste0(project_location, "/summary_analysis"))) {
dir.create(paste0(project_location, "/summary_analysis"))
}
set.seed(3221984)df_agg <- read.csv(paste0(project_location, '/summary_output/aggregate_data.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
clif_hospital_description <- read.csv(paste0(project_location, '/clif_hospital_description.csv'))
#Merge Hospital Data Into DF AGG
df_agg <- df_agg |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Do this Before Filtering
df_unique_hosp <- df_agg |> distinct(hospital_id, Hospital_Type)
cat('\nHow Many Academic v Community Hospitals are In THese Data')##
## How Many Academic v Community Hospitals are In THese Datatable(df_unique_hosp$Hospital_Type)##
## Academic Community
## 12 25df_hosptype_summary <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.df_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 13 10 12
## 2 Post-COVID 12 11 11
## 3 Pre-COVID 8 8 10#Filter Out Hospitals with Less Than 10 Patients
df_agg <- df_agg |>
#Filter OUT HOspitals-Periods with < 10 observations
filter(n_patients>=10)
#Now After Filtering
cat('\nAfer Filtering Out Hospitals with >= 10 Patients \nHow Many Academic v Community Hospitals are In THese Data')##
## Afer Filtering Out Hospitals with >= 10 Patients
## How Many Academic v Community Hospitals are In THese Datatable(df_agg$Hospital_Type)##
## Academic Community
## 33 50df_hosptype_summary_post <- df_agg |>
group_by(Hospital_Type, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, Hospital_Type)## `summarise()` has grouped output by 'Hospital_Type'. You can override using the
## `.groups` argument.cat('\nProning Counts in Hospitals with >= 10 Patients')##
## Proning Counts in Hospitals with >= 10 Patientsdf_hosptype_summary## # A tibble: 6 × 6
## # Groups: Hospital_Type [2]
## Hospital_Type study_period n_Hospitals n_patients n_proned percent_proned
## <chr> <chr> <int> <int> <int> <dbl>
## 1 Academic COVID 12 1884 680 36.1
## 2 Community COVID 23 1085 643 59.3
## 3 Academic Post-COVID 11 1098 168 15.3
## 4 Community Post-COVID 23 498 150 30.1
## 5 Academic Pre-COVID 10 982 66 6.72
## 6 Community Pre-COVID 16 213 30 14.1#Hospital Size by Period
df_hospsize <- df_agg |>
group_by(study_period) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description')##
## Hospital Size Descriptiondf_hospsize## # A tibble: 3 × 4
## study_period small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 COVID 12 10 12
## 2 Post-COVID 7 11 11
## 3 Pre-COVID 3 7 10#Academic Size Summary
df_hospsize <- df_agg |>
group_by(Hospital_Type) |>
summarise(
small_icus =sum(icu_bed_cat=='small'),
medium_icus = sum(icu_bed_cat=='medium'),
large_icus = sum(icu_bed_cat=='large')
)
cat('\nHospital Size Description by Hospital Type')##
## Hospital Size Description by Hospital Typedf_hospsize## # A tibble: 2 × 4
## Hospital_Type small_icus medium_icus large_icus
## <chr> <int> <int> <int>
## 1 Academic 0 3 30
## 2 Community 22 25 3#Proning By Hospital Size
df_hospsize <- df_agg |>
group_by(icu_bed_cat, study_period) |>
summarise(
n_Hospitals=n(),
n_patients=sum(n_patients, na.rm=T),
n_proned=sum(observed_prone, na.rm=T),
percent_proned=round((n_proned/n_patients)*100, 2)
) |>
arrange(study_period, icu_bed_cat)## `summarise()` has grouped output by 'icu_bed_cat'. You can override using the
## `.groups` argument.cat('\nProning Counts by Hospital Size and Period')##
## Proning Counts by Hospital Size and Perioddf_hospsize## # A tibble: 9 × 6
## # Groups: icu_bed_cat [3]
## icu_bed_cat study_period n_Hospitals n_patients n_proned percent_proned
## <fct> <chr> <int> <int> <int> <dbl>
## 1 small COVID 12 433 300 69.3
## 2 medium COVID 10 623 283 45.4
## 3 large COVID 12 1906 738 38.7
## 4 small Post-COVID 7 128 53 41.4
## 5 medium Post-COVID 11 335 77 23.0
## 6 large Post-COVID 11 1109 179 16.1
## 7 small Pre-COVID 3 43 6 14.0
## 8 medium Pre-COVID 7 156 23 14.7
## 9 large Pre-COVID 10 979 67 6.84#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned <- df_agg
df_agg_proned <- df_agg_unproned %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned <- df_agg_unproned |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg <- bind_rows(df_agg_proned, df_agg_unproned)
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
global_agg_byoutcome_period <- read.csv(paste0(project_location, '/summary_output/global_aggregate_byoutcome_period.csv')) |>
mutate(prone_log_odds = log(prone_rate_adjust / (1 - prone_rate_adjust))) |>
select(hospital_id, prone_12hour_outcome, prone_log_odds, study_period) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Back with DF Agg
df_agg <- left_join(df_agg, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, study_period) |>
arrange(hospital_id, study_period, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, study_period, prone_12hour_outcome)`#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Also Includes Hospital Type and Hospital Type and Period Interaction
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_full <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_full)
summary(agg_brm_period_full)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + study_period * Hospital_Type + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.91 0.13 0.70
## sd(study_periodPostMCOVID) 1.01 0.18 0.72
## sd(study_periodPreMCOVID) 1.00 0.27 0.58
## cor(study_periodCOVID,study_periodPostMCOVID) 0.77 0.11 0.51
## cor(study_periodCOVID,study_periodPreMCOVID) 0.59 0.21 0.11
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.52 0.23 0.03
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.18 1.00 1198
## sd(study_periodPostMCOVID) 1.41 1.00 1493
## sd(study_periodPreMCOVID) 1.61 1.00 1583
## cor(study_periodCOVID,study_periodPostMCOVID) 0.93 1.00 1605
## cor(study_periodCOVID,study_periodPreMCOVID) 0.92 1.00 1599
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.90 1.00 1461
## Tail_ESS
## sd(study_periodCOVID) 2125
## sd(study_periodPostMCOVID) 2311
## sd(study_periodPreMCOVID) 2425
## cor(study_periodCOVID,study_periodPostMCOVID) 2070
## cor(study_periodCOVID,study_periodPreMCOVID) 1412
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1962
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI
## study_periodCOVID -0.11 0.26 -0.61
## study_periodPostMCOVID -1.19 0.31 -1.82
## study_periodPreMCOVID -2.11 0.34 -2.80
## Hospital_TypeCommunity 0.85 0.32 0.20
## study_periodPostMCOVID:Hospital_TypeCommunity 0.10 0.31 -0.52
## study_periodPreMCOVID:Hospital_TypeCommunity -0.16 0.49 -1.11
## u-95% CI Rhat Bulk_ESS Tail_ESS
## study_periodCOVID 0.40 1.01 503 873
## study_periodPostMCOVID -0.60 1.01 716 1442
## study_periodPreMCOVID -1.48 1.00 1242 2148
## Hospital_TypeCommunity 1.48 1.01 552 1013
## study_periodPostMCOVID:Hospital_TypeCommunity 0.74 1.00 1872 2212
## study_periodPreMCOVID:Hospital_TypeCommunity 0.82 1.00 1797 1795
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## 3.26 seconds (Total)
## Chain 4:
## Chain 2: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 2.033 seconds (Warm-up)
## Chain 3: 2.35 seconds (Sampling)
## Chain 3: 4.383 seconds (Total)
## Chain 3:
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 2.847 seconds (Warm-up)
## Chain 2: 2.309 seconds (Sampling)
## Chain 2: 5.156 seconds (Total)
## Chain 2:#Now Generate Odds Ratios for Communtiy vs Academic Proning by Period
full_posterior <- as_draws_df(agg_brm_period_full) |>
mutate(pre_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPreMCOVID:Hospital_TypeCommunity`),
covid = exp(b_Hospital_TypeCommunity),
post_covid = exp(b_Hospital_TypeCommunity + `b_study_periodPostMCOVID:Hospital_TypeCommunity`))
summary_posterior <- data.frame(
variable = c(
'Pre-COVID Community Hospital',
'COVID Community Hospital',
'Post-COVID Community Hospital'
),
median_or = c(median(full_posterior$pre_covid),
median(full_posterior$covid),
median(full_posterior$post_covid)),
conf.low = c(quantile(full_posterior$pre_covid, p=0.025),
quantile(full_posterior$covid, p=0.025),
quantile(full_posterior$post_covid, p=0.025)),
conf.high = c(quantile(full_posterior$pre_covid, p=0.975),
quantile(full_posterior$covid, p=0.975),
quantile(full_posterior$post_covid, p=0.975)))
print(summary_posterior)
write.csv(summary_posterior, paste0(project_location, '/summary_output/summary_posterior.csv'))####Predictions Command from Marginal Effects with 4000 posterior draws
reliability_adjust_full <- predictions(agg_brm_period_full,
type = 'response', ndraws=4000, re_formula = ~ (0 + study_period | hospital_id))
reliability_adjust_full <- as_tibble(reliability_adjust_full) |>
left_join(df_agg %>%
select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds', 'prone_log_odds', 'study_period', 'prone_12hour_outcome'),
by = c('hospital_id', 'study_period', 'prone_log_odds')) |>
#Need to Recombine Prone_12hour_outcome Estimates into 1 Per Hospital and Period
group_by(hospital_id, study_period) |>
mutate(
estimate_add=sum(estimate, na.rm=T),
ci_low_add=sum(conf.low, na.rm=T),
ci_hi_add=sum(conf.high, na.rm=T),
n_patients=sum(n_patients, na.rm=T)
) |>
arrange(hospital_id, study_period, -prone_12hour_outcome) |> #THis Arranges so Row Number 1 is the prone_12hour_outcome==1
filter(row_number()==1) |>
ungroup() |>
mutate(adjust_rate=estimate_add/n_patients) |>
mutate(ci_low=ci_low_add/n_patients) |>
mutate(ci_hi=ci_hi_add/n_patients) |>
mutate(period_rank=fcase(
study_period=='Pre-COVID', 0,
study_period=='COVID', 1,
study_period=='Post-COVID', 2
)) |>
#Orders for Graph
arrange(period_rank, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
mutate(actual_rate=observed_prone/n_patients) |>
mutate(grouping='Actual Rate') |>
left_join(df_agg %>% select('hospital_id', 'Hospital_Type', 'Number_of_ICU_Beds')) |>
distinct()## Warning in left_join(as_tibble(reliability_adjust_full), df_agg %>% select("hospital_id", : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 13 of `x` matches multiple rows in `y`.
## ℹ Row 13 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.## Joining with `by = join_by(hospital_id, Number_of_ICU_Beds)`## Warning in left_join(mutate(mutate(mutate(arrange(mutate(mutate(mutate(mutate(ungroup(filter(arrange(mutate(group_by(left_join(as_tibble(reliability_adjust_full), : Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 1 of `x` matches multiple rows in `y`.
## ℹ Row 151 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.break_points <- reliability_adjust_full |>
mutate(pre_covid_num=fcase(
row_number()!=1 & study_period=='Pre-COVID' & lead(study_period)=='COVID', hospital_rank),
post_covid_num=fcase(
study_period=='COVID' & lead(study_period)=='Post-COVID', hospital_rank
)) |>
summarise(
pre_covid_num=mean(pre_covid_num,na.rm=T)+0.5,
post_covid_num=mean(post_covid_num, na.rm=T)+0.5)
caterpillar_hospital_type <- ggplot(reliability_adjust_full,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.1),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals",
y = "Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
geom_vline(xintercept = break_points$pre_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 7, y = 0.95, label = 'Pre Pandemic', fontface = 'bold') +
geom_vline(xintercept = break_points$post_covid_num, linetype = 2, linewidth=0.3) +
annotate("text", x = 35.5, y = 0.95, label = 'Pandemic', fontface = 'bold') +
annotate("text", x = 68, y = 0.95, label = 'Post Pandemic', fontface = 'bold') +
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.5))
caterpillar_hospital_type
ggsave('caterpillar_hospital_byperiod_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
#With Actual Rate
#Add the Actual Value
with_rate <- caterpillar_hospital_type +
geom_point(aes(y=observed_prone/n_patients), color = 'maroon', alpha = 0.25)
with_rate
ggsave('caterpillar_fully_adjusted.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)#Stack Caterpillar Plots Vertically
vertical_data <- reliability_adjust_full |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=row_number()) |>
ungroup()
vertical <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.1, 0.90),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 9)) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1,
scales='free_x')## Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
## 3.5.0.
## ℹ Please use the `legend.position.inside` argument of `theme()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.ggsave('vertical_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
verticald <-ggplot(vertical_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period), color='lightgrey') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type)) +
geom_point(aes(color=Hospital_Type)) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
axis.line.x = element_line(color = "black", linewidth = 0.25),
strip.background = element_rect(fill = "grey90", color = NA),
strip.text = element_text(color = "black"),
legend.position = c(0.75, 0.85),
legend.background = element_rect(fill = "white", color = "grey"),
legend.title = element_text(size = 12),
legend.key.height = unit(1.5, "lines"),
legend.key.width = unit(2.2, 'lines')) +
facet_wrap(~factor(study_period, c('Pre-COVID', 'COVID', 'Post-COVID')),
ncol=1)
ggsave('vertical_caterpillar_option_D.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)
stacked_data <- vertical_data |>
group_by(study_period) |>
arrange(study_period, adjust_rate) |>
mutate(hospital_rank=if_else(
study_period=='COVID', row_number(), NA_real_
)) |>
group_by(hospital_id) |>
fill(hospital_rank, .direction = 'updown')
stacked <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, color = study_period)) +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= study_period), alpha = 0.4) +
geom_point(aes(color=study_period, shape=Hospital_Type), size = 2.5) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Hokusai3") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank())
ggsave('stacked_caterpillar.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked2 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, linetype=study_period), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) +
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Demuth") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by COVID Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals"
) +
theme_minimal() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar2.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).stacked_data <- stacked_data |>
group_by(study_period) |>
mutate(period_n=sum(n_patients)) |>
ungroup()
stacked3 <- ggplot(stacked_data,
aes(x = hospital_rank,
y = adjust_rate)) +
geom_line(aes(x = hospital_rank,
y = adjust_rate, group=study_period, alpha=period_n), color = 'black') +
geom_linerange(aes(ymin = ci_low, ymax = ci_hi, color= Hospital_Type), linewidth=0.75, alpha = 0.4) + scale_alpha(guide="none")+
geom_point(aes(color=Hospital_Type, shape=study_period), size = 2.25) + # Add points for center effects
scale_x_continuous(
breaks = seq(1, nrow(reliability_adjust_full), by=1),
labels = NULL # Optionally use hospital IDs as labels
) +
scale_y_continuous(breaks=seq(0,1.0, by=0.2),
labels = scales::percent,
limits = c(0,1.0)) +
MetBrewer::scale_color_met_d("Egypt") + #MetBrewer has many pallettes to choose from
labs(
x = "Hospitals\nRanked by Pandemic Period Proning Use",
y = "% Eligible Patients Proned Within 12 Hours",
caption = "Error bars represent 95% credible intervals",
shape="Study Period",
color= "Hospital Type"
) +
scale_shape_manual(values= c("circle", "triangle", "square"), labels=c("Pandemic", "Post-pandemic", "Pre-pandemic"))+
theme_classic() +
theme(panel.grid.major.x=element_blank(),
panel.grid.minor.x=element_blank(),
panel.grid.minor.y=element_blank(),
legend.position = c(0.15, 0.75),
legend.background = element_rect(fill = "white", color = "grey"))
ggsave('stacked_caterpillar3.pdf',
path = paste0(project_location, '/summary_output/graphs/'),
width = 12, # adjust width to make it wide
height = 6, # adjust height for PowerPoint-friendly proportions
units = "in") # units set to inches for precise control)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_segment()`).## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).posterior <- as_draws_df(agg_brm_period_full)
# Extract the standard deviations (SDs) for each study_period random effect
sd_pre <- posterior$`sd_hospital_id__study_periodPreMCOVID`
sd_covid <- posterior$`sd_hospital_id__study_periodCOVID`
sd_post <- posterior$`sd_hospital_id__study_periodPostMCOVID`
#3 Calculate MOR for Each Period
# MOR calculation function
calc_mor <- function(sd_samples) {
exp(sqrt(2) * sd_samples * qnorm(0.75)) # qnorm(0.75) ~ 0.674
}
# Compute MORs
mor_pre <- calc_mor(sd_pre)
mor_covid <- calc_mor(sd_covid)
mor_post <- calc_mor(sd_post)
# Differences Using The Latter Periods as Reference
diff_covid_pre <- mor_covid - mor_pre
diff_post_pre <- mor_post - mor_pre
diff_post_covid <- mor_post - mor_covid
# Posterior probabilities that the Difference in mOR is > 0 (Bayesian 'p-values')
prop_covid_gt_pre <- mean(diff_covid_pre > 0)
prop_post_gt_pre <- mean(diff_post_pre > 0)
prop_post_gt_covid <- mean(diff_post_covid > 0)
# Credible intervals
cr_covid_pre <- quantile(diff_covid_pre, c(0.025, 0.5, 0.975))
cr_post_pre <- quantile(diff_post_pre, c(0.025, 0.5, 0.975))
cr_post_covid <- quantile(diff_post_covid, c(0.025, 0.5, 0.975))
# Summarize
mor_summary <- data.frame(
Period = c("Pre-MCOVID", "COVID", "Post-MCOVID"),
MOR_Mean = c(mean(mor_pre), mean(mor_covid), mean(mor_post)),
MOR_Median = c(median(mor_pre), median(mor_covid), median(mor_post)),
MOR_2.5 = c(quantile(mor_pre, 0.025), quantile(mor_covid, 0.025), quantile(mor_post, 0.025)),
MOR_97.5 = c(quantile(mor_pre, 0.975), quantile(mor_covid, 0.975), quantile(mor_post, 0.975))
)
print(mor_summary)## Period MOR_Mean MOR_Median MOR_2.5 MOR_97.5
## 1 Pre-MCOVID 2.690205 2.517063 1.745838 4.657248
## 2 COVID 2.396583 2.354657 1.944955 3.090758
## 3 Post-MCOVID 2.670977 2.587785 1.990022 3.836999#Summarize the Differences in Median Odds Ratios by Period
diff_mor <- data.frame(
contrast = c('COVID-Pre', 'Post-COVID', 'Post-Pre'),
MOR_diff_mean = c(mean(diff_covid_pre), mean(diff_post_covid), mean(diff_post_pre)),
MOR_diff_median = c(median(diff_covid_pre), median(diff_post_covid), median(diff_post_pre)),
MOR_diff_2.5 = c(quantile(diff_covid_pre, 0.025), quantile(diff_post_covid, 0.025), quantile(diff_post_pre, 0.025)),
MOR_diff_97.5 = c(quantile(diff_covid_pre, 0.975), quantile(diff_post_covid, 0.975), quantile(diff_post_pre, 0.975)),
prob_diff_gt_0 = c(prop_covid_gt_pre, prop_post_gt_covid, prop_post_gt_pre)
)
print(diff_mor)## contrast MOR_diff_mean MOR_diff_median MOR_diff_2.5 MOR_diff_97.5
## 1 COVID-Pre -0.29362135 -0.15911771 -2.2612521 0.8422212
## 2 Post-COVID 0.27439383 0.22563714 -0.5031843 1.3230072
## 3 Post-Pre -0.01922753 0.07103335 -2.0632630 1.4426286
## prob_diff_gt_0
## 1 0.39675
## 2 0.73075
## 3 0.55000#Visualize
mor_df <- data.frame(
MOR = c(mor_covid, mor_pre, mor_post),
Period = factor(rep(c("COVID", "Pre-MCOVID", "Post-MCOVID"), each = length(mor_covid)))
)
ggplot(mor_df, aes(x = MOR, fill = Period)) +
geom_density(alpha = 0.5) +
labs(
title = "Period-Specific Median Odds Ratios (MOR \n (Fully Adjusted))",
x = "Median Odds Ratio (MOR)",
y = "Density"
) +
scale_x_continuous(seq(1,20, by=1), limits =c(1,20),
name='Distribution of Median Odds Ratos by Study Period') +
theme_minimal()
ggsave('mor_distribution_fully.pdf',
path = paste0(project_location, '/summary_output/graphs/'))## Saving 7 x 5 in image#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#Adjusted for COVID Status
#First Prepare Dataset
#Open Aggregate Data with COVID
df_agg_covid <- read.csv(paste0(project_location, '/summary_output/aggregate_data_wcovid.csv')) |>
#Prepare to merge in Hospital Details
mutate(
hospital_id=fifelse(
hospital_id=='HUP', 'penn_hup', hospital_id),
hospital_id=fifelse(
hospital_id=='HUPME', 'penn_hupme', hospital_id),
hospital_id=fifelse(
hospital_id=='PAH', 'penn_pah', hospital_id),
hospital_id=fifelse(
hospital_id=='PPMC', 'penn_ppmc', hospital_id),
hospital_id=fifelse(
hospital_id=='UMCOP', 'penn_umcop', hospital_id),
hospital_id=fifelse(
hospital_id=='CCH', 'penn_cch', hospital_id)
)
#Merge Hospital Data Into DF AGG
df_agg_covid <- df_agg_covid |>
left_join(clif_hospital_description %>% select(hospital_id, Hospital_Type, Number_of_ICU_Beds)) |>
#Define Hospital_Size Categorical Variable
mutate(icu_bed_cat=fcase(
Number_of_ICU_Beds<20, 0,
Number_of_ICU_Beds>=20 & Number_of_ICU_Beds<50, 1,
Number_of_ICU_Beds>=50, 2
)) |>
## 1. numeric 0 / 1 / 2 bucket --------------------------------------------
mutate(
icu_bed_cat = fcase(
Number_of_ICU_Beds < 20, 0L, # small
Number_of_ICU_Beds >= 20 & Number_of_ICU_Beds < 50, 1L, # medium
Number_of_ICU_Beds >= 50, 2L, # large
default = NA_integer_
)
) %>%
## 2. labelled factor ------------------------------------------------------
mutate(
icu_bed_cat = factor(
icu_bed_cat,
levels = 0:2, # make sure levels match
labels = c("small", "medium", "large")
)
)## Joining with `by = join_by(hospital_id)`#Create Row For Each Proned Summary Data and Unproned Summary Data
df_agg_unproned_covid <- df_agg_covid
df_agg_proned_covid <- df_agg_unproned_covid %>%
mutate(n_patients = observed_prone, prone_12hour_outcome = 1)
df_agg_unproned_covid <- df_agg_unproned_covid |>
mutate(n_patients = not_prone, observed_prone = 0, prone_12hour_outcome = 0)
df_agg_covid <- bind_rows(df_agg_proned_covid, df_agg_unproned_covid) |>
#Only have the propensity at the period and outcome level so need to add period back here
mutate(study_period=fcase(
period_sarscov2=='COVID_SarsCov2Neg-Unk', 'COVID',
period_sarscov2=='COVID_SarsCov2Pos', 'COVID',
period_sarscov2=='Post-COVID_SarsCov2Neg-Unk', 'Post-COVID',
period_sarscov2=='Post-COVID_SarsCov2Pos', 'Post-COVID',
period_sarscov2=='Pre-COVID', 'Pre-COVID'
))
#Merge in Proning Adjusted Rate for Each Hospital By Proned/Unproned Status
#Merge Back with DF Agg COVID
df_agg_covid <- left_join(df_agg_covid, global_agg_byoutcome_period) |>
#IF there were no patients proned can use the popensity for unproned patients by filling here
group_by(hospital_id, period_sarscov2) |>
arrange(hospital_id, period_sarscov2, prone_12hour_outcome) |>
fill(prone_log_odds, .direction = 'down') |>
ungroup()## Joining with `by = join_by(hospital_id, prone_12hour_outcome, study_period)`set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PreMCOVID'),
#COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2Pos'),
prior(normal(0, 1), class = "b", coef = 'period_sarscov2COVID_SarsCov2NegMUnk'),
#Post-COVID Period Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2Pos'),
prior(normal(-2, 1), class = "b", coef = 'period_sarscov2PostMCOVID_SarsCov2NegMUnk'),
#Flat Priors for Hospital Type and Size
prior(normal(0,1), class = "b", coef = 'Hospital_TypeCommunity'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_covid <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type +
(0 + period_sarscov2 | hospital_id)),
data = df_agg_covid,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Warning: Rows containing NAs were excluded from the model.## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start samplingplot(agg_brm_period_covid)
summary(agg_brm_period_covid)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + period_sarscov2 + Hospital_Type + (0 + period_sarscov2 | hospital_id)
## Data: df_agg_covid (Number of observations: 318)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 37)
## Estimate
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.94
## sd(period_sarscov2COVID_SarsCov2Pos) 0.91
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.97
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.61
## sd(period_sarscov2PreMCOVID) 0.96
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.86
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.79
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.68
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.56
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.52
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.62
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.48
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.57
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.42
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.39
## Est.Error
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.15
## sd(period_sarscov2COVID_SarsCov2Pos) 0.12
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.16
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.28
## sd(period_sarscov2PreMCOVID) 0.25
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.09
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.12
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.13
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.27
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.26
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.26
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.21
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.19
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.22
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.30
## l-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 0.67
## sd(period_sarscov2COVID_SarsCov2Pos) 0.69
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.70
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 0.10
## sd(period_sarscov2PreMCOVID) 0.57
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.64
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.51
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.37
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.09
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) -0.10
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) -0.05
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.01
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.14
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) -0.05
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) -0.26
## u-95% CI
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.28
## sd(period_sarscov2COVID_SarsCov2Pos) 1.18
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.31
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.23
## sd(period_sarscov2PreMCOVID) 1.54
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 0.97
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.95
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 0.88
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.92
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 0.90
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 0.94
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.84
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.88
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 0.80
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 0.87
## Rhat
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2COVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## sd(period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 1.00
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1.00
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1.00
## Bulk_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 1395
## sd(period_sarscov2COVID_SarsCov2Pos) 1947
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2339
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 2033
## sd(period_sarscov2PreMCOVID) 1923
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1359
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 1448
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2034
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 3045
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 3958
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 3342
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 1887
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 2653
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2735
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 1198
## Tail_ESS
## sd(period_sarscov2COVID_SarsCov2NegMUnk) 2121
## sd(period_sarscov2COVID_SarsCov2Pos) 2828
## sd(period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2864
## sd(period_sarscov2PostMCOVID_SarsCov2Pos) 1504
## sd(period_sarscov2PreMCOVID) 2689
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2COVID_SarsCov2Pos) 1861
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 2282
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2NegMUnk) 3171
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2448
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PostMCOVID_SarsCov2Pos) 2573
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PostMCOVID_SarsCov2Pos) 2498
## cor(period_sarscov2COVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 2228
## cor(period_sarscov2COVID_SarsCov2Pos,period_sarscov2PreMCOVID) 2631
## cor(period_sarscov2PostMCOVID_SarsCov2NegMUnk,period_sarscov2PreMCOVID) 3061
## cor(period_sarscov2PostMCOVID_SarsCov2Pos,period_sarscov2PreMCOVID) 2439
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI
## period_sarscov2COVID_SarsCov2NegMUnk -0.99 0.22 -1.43 -0.56
## period_sarscov2COVID_SarsCov2Pos 0.46 0.22 0.01 0.88
## period_sarscov2PostMCOVID_SarsCov2NegMUnk -1.25 0.24 -1.71 -0.78
## period_sarscov2PostMCOVID_SarsCov2Pos -0.27 0.27 -0.83 0.25
## period_sarscov2PreMCOVID -2.29 0.29 -2.89 -1.78
## Hospital_TypeCommunity 0.84 0.25 0.34 1.33
## Rhat Bulk_ESS Tail_ESS
## period_sarscov2COVID_SarsCov2NegMUnk 1.01 753 1254
## period_sarscov2COVID_SarsCov2Pos 1.01 846 1578
## period_sarscov2PostMCOVID_SarsCov2NegMUnk 1.00 975 1800
## period_sarscov2PostMCOVID_SarsCov2Pos 1.00 1768 2496
## period_sarscov2PreMCOVID 1.00 1617 2477
## Hospital_TypeCommunity 1.01 973 1833
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).#Compare Differences in COVID vs Not COVID in COVID and POST COVID Period
post_covid <- as_draws_df(agg_brm_period_covid)
# Compute OR for the contrast
delta_covid_period <- post_covid$b_period_sarscov2COVID_SarsCov2Pos - post_covid$b_period_sarscov2COVID_SarsCov2NegMUnk
or <- exp(delta_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period \n', posterior_summary)## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During COVID Period
## 4.2533 3.20544 5.561699# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2Pos - post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for COVID vs Not COVID During Post-COVID Period
## 2.666524 1.570643 4.490486#COVID Negative in COVID Period vs Pre-COVID
# Compute OR for the contrast Of COVID vs NOT COVID in POST COVID Period
delta_post_covid_period <- post_covid$b_period_sarscov2PostMCOVID_SarsCov2NegMUnk - post_covid$b_period_sarscov2PreMCOVID
or <- exp(delta_post_covid_period)
# Summarize
posterior_summary <- quantile(or, probs = c(0.5, 0.025, 0.975))
cat('\nOdds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID \n', posterior_summary)##
## Odds Ratio and 95 Credible Interview for SARS-CoV2 Neg in COVID vs Pre-COVID
## 2.812252 1.649069 5.135836#ICU bed capacity Not Included as an additional covariable (in addition to hospital type becaues of co-linearity)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
#ICU Capacity is colinear with academic or community so have taken that out of the model
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% range
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Flat Priors for Hospital Size
prior(normal(0,1), class = "b", coef = 'icu_bed_catlarge'),
prior(normal(0,1), class = "b", coef = 'icu_bed_catmedium'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period and allowing for interaction term between hospital type and period
agg_brm_period_size <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + icu_bed_cat +
(0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling# Examine model fit and summaries
plot(agg_brm_period_size)
summary(agg_brm_period_size)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + icu_bed_cat + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.93 0.13 0.71
## sd(study_periodPostMCOVID) 1.03 0.18 0.72
## sd(study_periodPreMCOVID) 1.02 0.26 0.60
## cor(study_periodCOVID,study_periodPostMCOVID) 0.80 0.10 0.54
## cor(study_periodCOVID,study_periodPreMCOVID) 0.63 0.19 0.18
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.54 0.22 0.06
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.24 1.00 1224
## sd(study_periodPostMCOVID) 1.42 1.00 1581
## sd(study_periodPreMCOVID) 1.61 1.00 1858
## cor(study_periodCOVID,study_periodPostMCOVID) 0.94 1.00 1732
## cor(study_periodCOVID,study_periodPreMCOVID) 0.93 1.00 1928
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.91 1.00 1403
## Tail_ESS
## sd(study_periodCOVID) 1816
## sd(study_periodPostMCOVID) 2556
## sd(study_periodPreMCOVID) 2671
## cor(study_periodCOVID,study_periodPostMCOVID) 2488
## cor(study_periodCOVID,study_periodPreMCOVID) 2396
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 2180
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.67 0.25 0.16 1.16 1.00 1174
## study_periodPostMCOVID -0.35 0.28 -0.90 0.19 1.00 1309
## study_periodPreMCOVID -1.42 0.33 -2.11 -0.78 1.00 1223
## icu_bed_catmedium -0.32 0.36 -1.00 0.39 1.00 1232
## icu_bed_catlarge -0.55 0.33 -1.19 0.08 1.00 977
## Tail_ESS
## study_periodCOVID 1917
## study_periodPostMCOVID 2112
## study_periodPreMCOVID 2316
## icu_bed_catmedium 2041
## icu_bed_catlarge 1667
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## ing)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 2.302 seconds (Warm-up)
## Chain 1: 2.729 seconds (Sampling)
## Chain 1: 5.031 seconds (Total)
## Chain 1:
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.822 seconds (Warm-up)
## Chain 4: 2.316 seconds (Sampling)
## Chain 4: 5.138 seconds (Total)
## Chain 4:
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 3.078 seconds (Warm-up)
## Chain 3: 3.011 seconds (Sampling)
## Chain 3: 6.089 seconds (Total)
## Chain 3:
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 3.266 seconds (Warm-up)
## Chain 2: 3.066 seconds (Sampling)
## Chain 2: 6.332 seconds (Total)
## Chain 2:#Now Generate Odds Ratios for Hospital Size Contrasts
full_posterior <- as_draws_df(agg_brm_period_size) |>
mutate(medium_v_small = exp(b_icu_bed_catmedium ),
large_v_small = exp(b_icu_bed_catlarge))
summary_posterior_size <- data.frame(
variable = c(
'Medium vs Small ICU Capacity',
'Large vs Small ICU Capacity'
),
median_or = c(median(full_posterior$medium_v_small),
median(full_posterior$large_v_small)),
conf.low = c(quantile(full_posterior$medium_v_small, p=0.025),
quantile(full_posterior$large_v_small, p=0.025)),
conf.high = c(quantile(full_posterior$medium_v_small, p=0.975),
quantile(full_posterior$large_v_small, p=0.975)))
print(summary_posterior_size)
write.csv(summary_posterior_size, paste0(project_location, '/summary_output/summary_posterior_size_analysis.csv'))#Use Methods from Yarnell et al (PMID: 31359459)
#Bayesian Hierarchial Model with Random Effects for Hospital and Study Period
set.seed(32284)
priors <- c(
#Pre-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% with 95% of prior probability in 5-50% rnage
prior(normal(-2, 1), class = "b", coef = 'study_periodPreMCOVID'),
#COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 50% in COVID era and with 95% of these probability ranging from 12-88%
prior(normal(0, 1), class = "b", coef = 'study_periodCOVID'),
#Post-COVID Proning Prior: Weakly Informative Centered Around Prob of Proning of 12% in COVID era and with 95% of these probability ranging from 5-50%
prior(normal(-2, 1), class = "b", coef = 'study_periodPostMCOVID'),
#Half Normal Weakly Informative but Regularizing Prior for Random Effects
prior(normal(0,1), class = "sd"))
#Bayesian Hierarchical Model with varying slopes by study_period
agg_brm_period_strat <- brm(
bf(observed_prone | trials(n_patients) ~
0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)),
data = df_agg,
family = binomial,
prior = priors,
cores = 4,
control = list(adapt_delta = 0.99)
)## Compiling Stan program...## Trying to compile a simple C file## Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
## using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.5)’
## using SDK: ‘’
## clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/Rcpp/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/unsupported" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/src/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppParallel/include/" -I"/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include '/Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp' -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o
## In file included from <built-in>:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
## In file included from /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
## /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
## 679 | #include <cmath>
## | ^~~~~~~
## 1 error generated.
## make: *** [foo.o] Error 1## Start sampling## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess# Examine model fit and summaries
plot(agg_brm_period_strat)
summary(agg_brm_period_strat)## Family: binomial
## Links: mu = logit
## Formula: observed_prone | trials(n_patients) ~ 0 + offset(prone_log_odds) + study_period + (0 + study_period | hospital_id)
## Data: df_agg (Number of observations: 166)
## Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup draws = 4000
##
## Multilevel Hyperparameters:
## ~hospital_id (Number of levels: 36)
## Estimate Est.Error l-95% CI
## sd(study_periodCOVID) 0.98 0.13 0.77
## sd(study_periodPostMCOVID) 1.09 0.18 0.79
## sd(study_periodPreMCOVID) 1.03 0.25 0.63
## cor(study_periodCOVID,study_periodPostMCOVID) 0.82 0.09 0.61
## cor(study_periodCOVID,study_periodPreMCOVID) 0.65 0.19 0.22
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.57 0.22 0.11
## u-95% CI Rhat Bulk_ESS
## sd(study_periodCOVID) 1.27 1.00 790
## sd(study_periodPostMCOVID) 1.50 1.00 1292
## sd(study_periodPreMCOVID) 1.59 1.00 1646
## cor(study_periodCOVID,study_periodPostMCOVID) 0.95 1.00 1480
## cor(study_periodCOVID,study_periodPreMCOVID) 0.94 1.00 1424
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 0.92 1.00 1306
## Tail_ESS
## sd(study_periodCOVID) 1734
## sd(study_periodPostMCOVID) 2526
## sd(study_periodPreMCOVID) 2829
## cor(study_periodCOVID,study_periodPostMCOVID) 1909
## cor(study_periodCOVID,study_periodPreMCOVID) 1739
## cor(study_periodPostMCOVID,study_periodPreMCOVID) 1794
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## study_periodCOVID 0.39 0.17 0.05 0.72 1.02 362
## study_periodPostMCOVID -0.64 0.20 -1.06 -0.26 1.01 688
## study_periodPreMCOVID -1.72 0.26 -2.24 -1.23 1.01 1103
## Tail_ESS
## study_periodCOVID 1009
## study_periodPostMCOVID 1469
## study_periodPreMCOVID 2331
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
## ration: 1000 / 2000 [ 50%] (Warmup)
## Chain 4: Iteration: 1001 / 2000 [ 50%] (Sampling)
## Chain 1: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 2: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 4: Iteration: 1200 / 2000 [ 60%] (Sampling)
## Chain 3: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 1: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 2: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 1400 / 2000 [ 70%] (Sampling)
## Chain 3: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 1: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 1:
## Chain 1: Elapsed Time: 1.733 seconds (Warm-up)
## Chain 1: 1.942 seconds (Sampling)
## Chain 1: 3.675 seconds (Total)
## Chain 1:
## Chain 2: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 2:
## Chain 2: Elapsed Time: 1.826 seconds (Warm-up)
## Chain 2: 1.936 seconds (Sampling)
## Chain 2: 3.762 seconds (Total)
## Chain 2:
## Chain 4: Iteration: 1600 / 2000 [ 80%] (Sampling)
## Chain 3: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 3:
## Chain 3: Elapsed Time: 1.791 seconds (Warm-up)
## Chain 3: 2.196 seconds (Sampling)
## Chain 3: 3.987 seconds (Total)
## Chain 3:
## Chain 4: Iteration: 1800 / 2000 [ 90%] (Sampling)
## Chain 4: Iteration: 2000 / 2000 [100%] (Sampling)
## Chain 4:
## Chain 4: Elapsed Time: 2.618 seconds (Warm-up)
## Chain 4: 1.928 seconds (Sampling)
## Chain 4: 4.546 seconds (Total)
## Chain 4:#Posterior Draws
sensitivity_posterior <- as_draws_df(agg_brm_period_strat) |>
mutate(covid_v_pre_covid=exp(b_study_periodCOVID-b_study_periodPreMCOVID),
covid_v_post_covid=exp(b_study_periodCOVID-b_study_periodPostMCOVID),
pre_covid_v_covid=exp(b_study_periodPreMCOVID-b_study_periodCOVID),
post_covid_v_covid=exp(b_study_periodPostMCOVID-b_study_periodCOVID),
post_covid_v_precovid=exp(b_study_periodPostMCOVID-b_study_periodPreMCOVID))
sensitivity_table <- data.frame(
variable = c('COVID_v_Pre-COVID',
'COVID_v_Post-COVID',
'Pre-COVID_v_COVID',
'Post-COVID_v_COVID',
'Post-COVID_v_Pre-COVID'),
or = c(median(sensitivity_posterior$covid_v_pre_covid),
median(sensitivity_posterior$covid_v_post_covid),
median(sensitivity_posterior$pre_covid_v_covid),
median(sensitivity_posterior$post_covid_v_covid),
median(sensitivity_posterior$post_covid_v_precovid)),
conf.low=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.025),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.025),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.025),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.025),
quantile(sensitivity_posterior$post_covid_v_precovid, p=0.025)),
conf.high=c(quantile(sensitivity_posterior$covid_v_pre_covid, p=0.975),
quantile(sensitivity_posterior$covid_v_post_covid, p=0.975),
quantile(sensitivity_posterior$pre_covid_v_covid, p=0.975),
quantile(sensitivity_posterior$post_covid_v_covid, p=0.975),
quantile(sensitivity_posterior$post_covid_v_precovid, p=0.975)))
sensitivity_table
write.csv(sensitivity_table, paste0(project_location, '/summary_output/sensitivity_table.csv'))